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ABSTRACT 

An intuitive trial-and-error approach to history-matching can be 
costly and time-consuming. Considerable efforts have therefore been 
made to automate history-matching procedures for implementation on high- 
Speed computers. 

Several such automatic history-matching methods have been reported, 
but, because their effectiveness and reliability remain to be 
demonstrated, they have not been well-received by practising reservoir 
engineers. 

This study examines presently available automatic history-matching 
algorithms, discusses those that appear most viable, and explores the 
fundamental weaknesses of these methods. In addition, the existence of 
alternative approaches to solving the inverse problem for analogous, but 


less complex, situations is noted. 
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NOMENCLATURE 


ENGLISH 

B = formation volume factor, m3/sm° 

qe es rock compressibility, Pant 

D = depth, measured vertically positive downwards, m 

g = acceleration due to gravity, m/s? 

k = absolute permeability, me 

Me = relative permeability, fraction 

Pego = gas-oil capillary pressure, Pa 

eo wineO nl -wacen capillary pressure, Pa 

Pr = phase pressure, Pa 

dp = phase fluid production(-)/injection(+) rate, sm?/s 
Re, = solubility of gas in liquid phase, sm?/sm° 

hae = solubility of gas in water phase, sm3/sm> 

So phase saturation, fraction 

t = time, s_- 

Vy = bulk volume, m2 

Kes = mole fraction of oil, water, or gas in aqueous phase 
eas mole fraction of oil, water, or gas in liquid phase 


a mole fraction of oil, water, or gaS in vapour phase 
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Greek 
At = time step, t™l- +t") 


Ax, Ay,AZ = block dimensions in x,y,z directions, m 


» = viscosity, Paes 

o = density, kg/m? 

= molar density, kg-mole/m° 
S = potential, Pa 

d = porosity, fraction 
Subscripts 

a = aqueous phase 

f = phase 

g = gas 


T,J.k = X;¥52- airection nodal subscripts 


1 = liquid phase 


0 = O71] 
v = vapour 
w = water 


MeVaz = %¥.2 GIbheceions 


Superscripts 


n = present time level index 


n+1 advanced time level index 
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CHAPTER ONE 
INTRODUCTION 

Increasingly sophisticated high-speed computing capabilities and 
attendant development of efficient numerical algorithms have made it 
possible to utilize complex mathematical models for analysis, 
development, and management of petroleum reservoirs. Such models can 
represent the simultaneous flow of 011, water, and gas, and the mutual 
interactions of these fluids, within the porous and permeable strata 
which comprise a reservoir. 

A reservoir simulator is used, first and foremost, for predicting 
the volume of oil recoverable by various production schemes. A wide 
spectrum of different operating parameters can thus be tested; and 
provided that the geology of the reservoir is known in sufficient 
detail, and a good estimate of the porosity and permeability of the rock 
can be made, the simulation results can then be used to formulate an 
optimal development strategy for a particular field. 

Rough estimates of the rock properties may be obtained from 
analysis of core samples taken from the reservoir. However, due to high 
costs of drilling on closely spaced centres - e.g., 15-20 m apart - the 
number of samples, is usually small; and for a large reservoir which is 
likely to be heterogeneous, such estimates are, at best, little more 
than educated guesses. Some refinement of the estimates is possible if 
other known data - i.e. the performance histories of the reservoir - are 
incorporated. Such so-called 'history-matching' (or ‘inverse 


simulation') involves using a reservoir simulator to ‘adjust’, within 
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Set bounds, future production/injection and pressure parameters to 
values consistent with previous reservoir behaviour. More specifically, 
it requires iterative computer simulations, with each run adjusting the 
primary input data (porosity, permeability, and relative permeability) 
until an acceptable 'match' has been abtained. 

While the concept of history-matching is straightforward, the 
achievement of a good match is difficult, primarily because of the non- 
uniqueness of the solution (see, e.g., Chapter 2). In practice, any 
history-match therefore seeks to determine the ‘best' set of input 
parameters within pre-determined bounds (which are usually based on 
estimates from core samples). 

Although history-matching is frequently performed by reservoir 
engineers on an intuitive, trial-and-error basis, considerable efforts 
have been made to design ‘automatic’ history-matching algorithms, i.e., 
to program, in some fashion, a computer to determine an optimal set of 
reservoir parameters. But despite these attempts to simplify (and speed 
up) the entire history-matching procedure, extant techniques appear not 
to have been well-received by practising reservoir engineers - possibly 
due to perceived shortcomings of such techniques. This study therefore 
examines the various, presently available, automatic history-matching 
algorithms, and evaluates their relative merits and weaknesses as viable 


alternatives to non-automatic techniques. 
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~ CHAPTER TWO 
REVIEW OF THE LITERATURE 

First attempts to estimate the geological properties of a petroleum 
reservoir for use in reservoir simulation studies were reported by 
Kruger (1961). He emphasized the importance of obtaining agreement 
between calculated and observed pressure histories, and suggested that 
reservoir properties be adjusted by small, successive perturbations, 
with a simulator run undertaken after each iteration in order to 
determine systematically an ‘acceptable’ match. 

Nelson (1962) provided a short analysis of the paper by Kruger 
(1961), and in particular, discussed conditions necessary to insure a 
unique and realistic determination of the permeability distribution. 

Efforts to automate this procedure, and thus exploit the speed and 
efficiency of computers, were initiated by Jacquard (1964), and were 
predicated on the mathematics of an electrical resistance - capacitance 
network (known as an ‘electric analyzer'). After demonstrating the 
close similarity between such a network and a reservoir model, Jacquard 
developed a method for interpreting pressure history data in terms of 
geological properties throughout a reservoir. The paper provides sample 
calculations of permeabilities for a simple, single-well reservoir with 
radial-circular heterogeneity, and indicates a pronounced sensitivity of 
results to minor variations in reservoir pressure. 

Jacquard & Jain (1965) extended this technique, with the aid of 
least-squares methods, to the case of a two-dimensional reservoir. For 


this purpose, the reservoir was divided into zones, each of which was 
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assumed to have a constant permeability. (Similar zonation was 
subsequently adopted by other authors of automatic history-matching 
algorithms in order to reduce the number of unknown parameters.) The 
study used a reservoir simulator to generate an ‘observed' pressure 
history, and then employed these data to compute reservoir parameters. 
It was found that, while retaining the original zonation led to 
excellent reconstitution of the reservoir, alternative zonations yielded 
poor estimates of reservoir properties. 

Dupuy (1968) pursued the work of Jacquard & Jain (1965), and found 
instabilities reflected in the fact that calculated permeabilities 
frequently oscillated from one iteration to the next. He also noted 
that certain permeabilities were unrealistically high, and suggested the 
solution technique emplayed be coupled with imposition of some limiting 
values. 

Jahns (1966), following the recommendations of Kruger (1961), 
Jacquard (1964), and Jacquard & Jain (1965), proposed the use of 
nonlinear regression analysis as means for selecting adjustment factors 
for reservoir parameters. An adaptation of the method of steepest 
descent, the technique was based on an analysis of effects that small, 
successive perturbations in each zone exert on performance histories 
predicted by a simulator. While the method was designed to handle 
Single-phase flow, Jahns pointed out that it would also be suitable for 
multiphase flow if saturation changes were relatively small. Jahns was, 
incidentally, the first to consider the adjustment of porosity as well 
as permeability, and further recognized the non-uniqueness of his 


geological estimates. The paper applied the algorithm to two actual 
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reservoirs; and, in both cases, the average computed properties compared 
favorably with the, albeit sparse, available data. It should be noted, 
however, that there was a tendency to predict rather extreme values - 
according to Jahns, due to the short-term performance history used in 
his calculations. 

Nelson (1968) approached the problem of determining reservoir 
permeability by incorporating an energy dissipation analysis. He 
examined the equation for single-phase flow, and reduced the resulting 
first-order partial differential equation in the unknown permeability to 
a system of characteristic equations; and this was then used to obtain a 
differential expression for permeability as a function of the known 
pressure distribution. For steady flow, for which the expression could 
be integrated, direct calculation of the permeability distribution along 
Successive streamlines was possible. However, although this method 
produced excellent results for several single-phase test reservoirs, it 
is now of little practical value: the need for additional estimation of 
porosity, modelling of three-dimensional systems - and, most important, 
modelling of multiphase flow (where relative permeabilities must also be 
estimated) - leads to mathematical equations that are too complex to be 
solved by this method. 

The work of Coats et al. (1970) led to the use of random selection 
of reservoir properties in a set of computer simulations. Assuming 
linear relationships between matching errors and reservoir properties, 
and calculating the coefficients for this relation by least-squares 
methods, these investigators employed linear programming to minimize the 


error and thus estimate both the reservoir porosity and permeability. 
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6. 
In several test cases, including a study of a two-phase reservoir, this 
procedure led to acceptable results. But although the system is 
basically nonlinear (Jacquard, 1964; Dupuy, 1968), the error-parameter 
relationship can be fairly well approximated with a linear function, 
especially where reservoirs possess relatively high porosities and 
permeabilities (Slater & Durrer, 1971; Carter et al., 1974). The 
principal advantage of the method of Coats et al. appears to lie in its 
ability to guarantee a global minimum error, as opposed to the local 
minima obtained by use of nonlinear optimization techniques. 

Basing their work on the same general principles as Jacquard 
(1964), Dupuy (1968), and Coats et al. (1970), Slater & Durrer (1971) 
used error-weighted gradients and a linear programming formulation to 
Systematically reduce differences between observed and calculated 
performance histories. In these studies, interference relationships 
were determined from the effects on performance histories of 
successively altering reservoir parameters. The interference 
relationships were then weighted and used to adjust the reservoir 
parameters. Slater & Durrer also discussed the impracticality of 
achieving a perfect match, especially when the accuracy of performance 
history measurements is questionable, and suggested reasonable 
guidelines for assessing when an acceptable match has been obtained. 
Examples involving hypothetical, single-phase reservoirs indicated a 
strong tendency to prediction of extreme values, and studies of actual 
reservoirs yielded poor results. Slater & Durrer consequently concluded 


that "in true reservoirs, where there exist uncertainties in pressure 
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the 
measurements and the choice of a simulation model, a matching process is 
much more complicated." 

A similar approach was employed by Thomas et al. (1972), who based 
their method on the classical Gauss-Newton least-squares procedure and 
parameter constraints, and made provision for handling highly nonlinear 
cases. They worked with data provided for three sample history- 
matchings used by Jahns (1966) and Coats et al. (1970) in order to 
compare methods and results, and found that while their algorithm in 
most cases required fewer simulator runs, it yielded matches almost 
identical to previously determined ones. 

Boberg et al. (1973), using a variation of the algorithm by Thomas 
et al. (1972), were able to improve the stability of the numerical 
solution, and reported reasonable results for a complex Middle East 
oilfield. 

Veatch & Thomas (1971) developed a novel direct method for handling 
the history-matching problem by treating the finite difference analogues 
of the partial differential equations for multiphase flow as a system of 
linear equations in the unknown reservoir properties. The technique is 
applicable to multiphase, compressible flow in heterogeneous reservoirs, 
and capable of calculating reservoir parameters in a single computer 
run. In addition, it contains provisions for constructing ‘best’ 
estimates where no unique solutions can be computed. However, 
relatively complete performance histories over varying periods of time 
are required as input data; and in most practical cases, where such data 
are lacking, it becomes necessary to employ interpolation schemes. 


Results of studies on several hypothetical reservoirs presented by 
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Veatch & Thomas were promising, but the method has not received 
widespread attention, and little commentary on the paper exists. 

Carter et al. (1974), primarily expanding the work of Jacquard 
(1964) and Jacquard & Jain (1965), focussed attention on the method by 
which differences between observed and calculated performance histories 
are used to adjust reservoir properties, and also introduced two 
modified, iterative linear programming procedures. As a result, a 
method was designed for handling history-matching in situations where 
more reservoir parameters must be determined than are observed (known as 
an ‘underdetermined' problem). Several examples of single-phase 
reservoirs were presented, and these indicated that linear programming 
formulations are able to yield reasonable estimates of reservoir 
parameters. For the case of underdetermined performance matching 
problem, a modified linear programming technique also produced 
acceptable results. 

Chavent et al. (1973) and Chen et al. (1974) approached automatic 
history-matching by making use of optimal control theory, and estimated 
reservoir properties by continuous functions rather than as discrete 
values. Case studies in both papers were limited to single-phase 
reservoirs, for which the necessary governing equations could be easily 
derived. The method compared favourably with standard, constant-zone 
gradient methods (such as those used by Slater & Durrer, 1971; Thomas et 
al., 1972; and Carter et al., 1974), and on average, reduced the 
required computer time. 

Wasserman et al. (1975) applied the optimal control methods of 
Chavent et al. (1973) and Chen et al. (1974) in an attempt to treat 


multiphase flow problems. Multiplication of the porosity and 
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permeability terms by saturation-dependent terms (obtained from running 
a multiphase simulator) seemed to transform, in effect, a multiphase 
reservoir to a 'pseudo' single-phase reservoir. However, while 
application of this technique to field reservoirs yielded acceptable 
results, none of the three studies cited above made any attempt to 
investigate the non-unique aspects of the solution in the context of 
optimal control theory. Nor did they discuss the likelihood of 
computing only a local minimum matching error. 

An initial study of history-matching in two-phase petroleum 
reservoirs was reported by van den Bosch & Seinfeld (1977), who 
investigated the feasibility of estimating two-phase reservoir 
properties, and considered a simple, hypothetical reservoir with radial- 
circular symmetry, incompressible flow, and a central producing well. 
Qualitative and quantitative results suggested that the ability to 
estimate reservoir parameters depends on the type of flow within the 
reservoir, and evidenced non-uniqueness in the solutions of the history- 
matching problem. 

A technique for history-matching based on estimation by Bayesian 
statistical methods was introduced by Gavalas et al. (1976). Analysis 
of a hypothetical, one-dimensional, single-phase reservoir indicated 
that estimation by such techniques produces better results than history- 
matching algorithms using zonation procedures. However, the accuracy of 
Bayesian estimates depends on the validity of the prior statistics (of 
reservoir properties) employed; and in practical situations, where 
detailed geological data are lacking, the procedure appears to be no 
more reliable, and probably less so, than earlier automatic history- 


matching methods. 
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Shah et al. (1978) presented a general analysis of errors arising 
in history-matching, with particular emphasis on error variation between 
Bayesian and constant-zone estimation techniques. From consideration of 
a hypothetical, one-dimensional, single-phase reservoir, they concluded 
that zonation techniques are clearly preferable when the location of 
zone boundaries is indicated by geological data. In cases where only 
Sparse data on reservoir properties are available, neither Bayesian nor 
zonation procedures can be employed with much confidence. 

Watson et al. (1980) further extended the optimal control algorithm 
of Chavent et al. (1973) and Chen et al. (1974) for automatic history- 
matching in two-phase reservoirs. They developed an algorithm for 
estimating relative permeability, as well as porosity and absolute 
permeability. Studies on hypothetical reservoirs produced reasonable 
estimates of reservoir properties. But the equations required for two- 
phase reservoir history-matching in this manner are extremely complex; 
and while it is theoretically possible to extend the algorithm to three- 
phase flow, practical considerations would seem to negate the 


feasibility of such a study. 
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CHAPTER THREE 
DEVELOPMENT OF A PETROLEUM RESERVOIR SIMULATOR 

Since a petroleum reservoir simulator is used as an essential 
component of automatic history-matching algorithms as well as for 
generating data that test the reliability of such algorithms, it is 
necessary to develop a suitable simulator before considering extant 
automatic history-matching techniques. This section discusses the 
construction of a three-phase (oil, water, gas), three-dimensional, 
compressible flow reservoir simulator in which oil is assumed to be non- 


volatible (such a simulator is known as a ‘black oil' simulator). 


Formulation of the Mathematical Model 

A mathematical model that simulates multiphase, multi-dimensional 
fluid flow in a petroleum reservoir is conveniently based on Darcy's 
equation of flow and on the law of conservation of mass (known as the 
‘equation of continuity') as applied to each phase. These governing 
relations have been established by Muskat (1949), Collins (1961), 
Scheidegger (1974), Thomas (1982), and others. 

Proceeding from these bases, a reservoir simulator for the present 
study was developed by assuming that 

(1) flow is laminar, viscous, and irrotational; 

(2) the flow process is isothermal ; 

(3) a thermodynamic equilibrium exists among the three phases; 

(4) relative permeability curves adequately depict multiphase 


flow; 


che 


abst 3 
rata yoha 
on? gsesucetb: ® 
Asfotahentbessitt mea | 
<n wt ol bemvege-ef Veneta ‘et notet umbe wae watt 
Jeaccoiuml "(ocsberd? 2 ak Puede zh vaverimee 


feralynpet tft lier ,eeenat sii eszalumie Peits *etj. Tabor teateenaiian'al oe 
pesabd. no: aga ivnalhavnds BPN WacSy iaverOndogia at al ate 
ar 2h Rwond) ‘ery %y nbbtavasener To Wet she nd bine wort: fe netneape « 
grinreci seed? yseerty MbeS.ot bel i gge 28 ('ytruni sds Fo cmpRueee 
EROL) wate fod: (REN) pedetMvo bone (dates nerd evel anole OS 
-enoidg one. (4801) comodT (2°22) senpebigeet: a 

gheseta eto? tosolumte sievis]es & p2szad seeds ton! paeheeoeh - 
jadi pntnuses yd begrisyeh 2ew ybuse 


piawatgesqisi Wh jpeitlosely ,1enimes 2h wott 

spautssityen | oe! epeso+g woly) She 

pReende swints of? gnome edetas mitrdi| uns aremoghomaty © 
sesriqghzium gofeed yistsupsbs eeveua vii (tésemeqg gvisbion 
swalt 


= 15 = 


ee 
(5) phase interchanges are restricted to gas solution and 
release from oil and water; and 
(6) the reservoir consists of a continuous porous medium which 
is externally bounded by an impermeable surface ( i.e., 
there is no fluid flow across the outer boundaries of the 
reservoir). 
The resulting simulator allows for effects of gravity, capillarity, 
fluid viscosity, relative permeabilities, gas solubility, and reservoir 
heterogeneity. 

Consider a three-dimensional reservoir of variable thickness in the 
x,y,z framework of a Cartesian coordinate system, and let this reservoir 
contain a distribution of line sources and sinks (i.e., production and 
injection wells, respectively) oriented parallel to the z-axis, whose 
strengths are given as functions of time and location. 

The general molar balance equation for each component i flowing 
in the oil (liquid), water (aqueous), and gas (vapour) phases through a 


porous medium is given by 
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Assuming that phase interchanges are restricted to the solution and 
release of gas from oil and water, the mole fractions Xoae Xove. Xwl 2 


and Xyy pare all zero. The non-zero mole fractions are given by 
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Substituting these values into the appropriate forms of Equation (3.1) 
yields the following phase equations for the simultaneous flow of oil, 


water, and gas: 
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Gas phase 
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If the porous medium of the reservoir is assumed to be fully 
Saturated by the three phases, and preferentially water wet at the oil- 


water interface, the phase saturations are related by 


S.a Sa he theca Ne (3.3) 
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The unknowns in the model are the pressure and saturations. Moreover, 
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containing a source or sink is specified, the other two can be 
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determined by making use of the relative mobility relationships 
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Thus, a model which describes the flow mechanics of a petroleum 
reservoir is therefore provided by using the coupled Equations (3.2) in 
combination with Equations (3.3) - (3.6) and appropriate initial, 
boundary, and source-sink conditions. (For further discussion on 
development of this basic model, see, for example, Peaceman, 1969; 


inomas.. 1982shor Bird etlal., 1960) 


Method of Solution 

Due largely to the nonlinear nature of the system of partial 
differential equations, no analytical solution of Equations (3.2) is 
possible (Jacquard, 1964; Dupuy, 1968). Consequently, the technique of 
finite differences is frequently employed to obtain a numerical 
solution. (Finite differences have been extensively used in the 
petroleum industry for solving problems in porous media flow, and have 
yielded results that stand in good accord with experimental data.) 

To utilize finite differences in the present case, the petroleum 
reservoir is defined by a block-centred grid that discretizes space 
variables, and the time continuum is subdivided into increments. Areal 
and cross-sectional sketches of a typical grid representation are shown 
in Figures 1 and 2. It is not necessary for the reservoir partitions to 
be equally spaced, and the blocks need not even be rectangular: 
depending upon the geological structure of the reservoir, circular, 
triangular, or curvilinear blocks may be used to more closely represent 


it (see Figure 3). 
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Figure 1. Rectangular Grid Representation of a Petroleum Reservoir 
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Figure 2. Cross-Sectional Grid Representation of a Petroleum Reservoir 
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Figure 3. 
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Cross-Sectional View. 
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Circular Grid Representation of a Petroleum Reservoir 
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The partial derivatives in Equations (3.2) are then expressed as 
algebraic finite difference equations for each grid-block; and after 
‘linearizing' them by freezing the coefficients over each time step, the 
resulting system of simultaneous equations can be solved repeatedly to 
obtain pressures, saturations, and flow rates at the desired advanced 
time levels and locations. Note that whereas an analytical solution 
requires the partial differential equations to be satisfied at every 
point in the domain, results of the finite difference approach are 
limited to pre-defined points in space and time, and subject to errors 
which depend on the magnitudes of the chosen increments, as well as 
Subject to truncation and roundoff errors incurred by the finite word 
length of the computer. 

The accuracy and effectiveness of the reservoir model is obviously 
improved if the number of grid-blocks is increased (and a better 
definition of the reservoir is provided). But a larger number of blocks 
rapidly increases computational difficulties; and given the heavy 
demands on machine time, it is essential to utilize the most efficient 
solution techniques. A single simulation may involve the solution of 
several hundred simultaneous equations up to several thousand times. 

Within the general framework of the finite difference approach, 
there exist two basic iterative methods of obtaining a numerical 
solution to the multiphase flow problem. The first of these, used in 
the present study, is the so-called ‘implicit pressure-explicit 
saturation’ method (hereafter referred to as the 'IMPES'). In this, all 
variables except one of the phase pressures (usually the oi] pressure) 


are eliminated by use of the phase, saturation and capillary pressure 
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equations. The resulting second-order equation, containing the unknown 
oil pressure, is then applied to each grid-block, and the system of 
equations thus obtained is solved for the oil pressure at the advanced 
time level for all blocks. Thereafter, rearrangement of the original 
partial differential equations for the oil and water phases yields the 
new oi] and water saturations at each block. The water and gas phase 
pressures are obtained from the capillary pressure equations, and gas 
Saturations are calculated from the saturation constraint equation. 

The other presently available solution method of which many 
variations exist, is known as the ‘implicit' or ‘simultaneous solution' 
method. This uses the saturation and capillary pressure equations to 
eliminate saturation terms from the right hand side of the three partial 
differential phase equations. The resulting three equations, which 
contain the phase pressures at the advanced time level as unknowns, are 
then written for each grid-block, and solved simultaneously. Note that 
while the implicit method is more stable (i.e., allows the use of a 
larger time increment) than the IMPES method, it requires a three times 
larger number of equations to be solved, and hence greatly increases the 
demand for computer time. (Chappelear & Rogers, 1973; Weinstein et al. 


1970). 


Derivation of the Numerical Model 


In order to translate the partial differential equations derived in 
Equation (3.2) to their algebraic finite difference analogues, it is 


convenient to rewrite them as follows: 
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The simplest method of ordering blocks in a three-dimensional grid, 
such as has been used in the present study, identifies each cell by the 
Subscripts 71, J]; and) k= (increasing im the normal “x-, y- and z- 
directions, respectively). Figures 4 and 5 illustrate this ordering 
scheme. 
Applying finite difference approximation techniques to Equations 


(3.7) then yields: 
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Figure 4. Conventional Identification of Blocks in Three-Dimensional 
Grid 
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Figure 5. Numbering System For Blocks in Three-Dimensional Grid. 
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where the difference operator (in any direction r, any block b) is 


defined by 


A (4.4.8) = M410 (Be 7-8) - 1/28, -Bp_7) 


Terms that account for the relative ease of fluid flow between grid- 


blocks are given by the 'transmissibilities'. For example 
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Also, for notational convenience, the symbols "'" and "g" are used to 


append transmissibility terms. For example, 


and 


To implement the IMPES solution method, the three phase equations 
(Equations (3.8)) are combined, with the aid of the saturation and 
capillary pressure equations (Equations (3.3) - (3.5)), in order to 
obtain a single equation in the unknown oil pressure. 

Consider, first, the expansion of the right hand side of the oil] 


phase equation (Equation (3.8a)): 
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Similarly, for the right hand side of the water phase equation (Equation 
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the right hand side of the gas phase equation (Equation (3.8c)) becomes 


So Sy 34 
Va lolRso B. i eu Bs eee) 


= 
a 


V \ 
b;,ontl ntl-,ntl,1l 1 , n+1 on 
Beat OG ? Le. E) aa ieplcrs suey) 
0 
n+1 
S n 
i Ke: ee kel oie n o ntl on 
ech ( o ) aaa 0 =3,) 
B B 
0 0 
nt1l-,nt1l,1 ,' ] n+1 on 
f ae hey teeta eas an llste Spay 


B 
W 


i, _ 
be : : 
bose! squtestelide Si 


aJ deme (a) ; 7 a9 
7 - : =) 


Sealy etew erie 26 wuts bner mibie ety vot erates: 


dotiaup®) nobteue 


-_ 


age See te ae ee Gepepre wk aay 
(ae) rl 


goats ,tne 


olay gh te ie 


genoosd ((986) notjsups) totfeups ‘seatq eg 947 TY) core ‘bred ftpin edd 


ght 1 é 
n —W Oran leet Nn o font] on 
wee gn ° ( oe e 0) aur an ow oon 
W W 
n+1 
S 
Nh nelvltyen sol gan g poten 
+ Sy 6 B) (Dreea a o (he en) 
g 
nel 
+ aS -s°)} ; (3.9c) 
B 
g 


Observe that all terms in these expansions are taken at any grid-block 
With andex(1,44K)< 

Substitution of Equations (3.4) and (3.5) into the above eliminates 
dependence on water and gas phase pressures at the n time level. To 
eliminate saturation terms at the n time level, Equations (3.8) 
combined with the expansions (3.9) are rearranged so that the right hand 
sides only contain terms in (sMtligny, The resulting right hand side 


expressions for each phase equation are then 
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Summing together with Equation (3.8c), and applying Equation (3.3) 


causes the right hand side terms to vanish: 
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Since calculating the changes in oi] pressures rather than the oil 
pressures themselves reduces round-off error, define 
n+1 n 
SPy ; See be. 4) aia 
15d ok Tusk 1,J5k 
and, in a notation convenient for computer implementation, the resulting 


sum of the three phase equations can then be written in the form 
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Writing Equation (3.10) for each grid-block, and ordering the grid 
as in Figure 6, then yields a system of linear equations of the form 
fixe =a. whene Ae 1Si the coett jerentamatrix containing 2.02. 0s jb, .r, 
H, S (see Figure 7), x is the column vector of the unknown OP 
terms, and b is the column vector containing the i,j,k terms. (For 
alternative ordering schemes, see Price & Coats, 1974). 

Once the phase pressures have been calculated for each grid-block, 
it is possible to calculate the phase saturations explicitly by means of 
the original difference equations. Rearrangement of the oil and water 


phase equations (Equations (3.8a) and (3.8b)) yields 
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Figure 6. Normal Ordering For Grid System 


Figure 7. Coefficient Matrix For Normally Ordered Three-Dimensional 
Grid System 
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The gas phase saturations can then be determined by means of Equation 
(3 63))" 

Since the coefficients of Equation (3.10) are dependent on both the 
unknown pressures and saturations, the solution process must be carried 
out by iteration. Initial estimates of pressures, saturations, and 
production/injection rates are made (usually based on graphical 
extrapolation schemes), and the resultant calculated values in each 
grid-block are compared with these estimates. If convergence to a 
solution (within a prescribed tolerance) has not been obtained, 
calculations are repeated with these latest values. Alternatively, a 
material balance is established, and if acceptable, the calculations 
advance to the next time step. If the material balance is unreasonable, 
the time step is diminished, and pressures, saturations, and 
production/injection rates are recalculated. The material balance 


check, defined as 


_ Net change in mass over time interval _ 1 


Me Net mass throughput 


») 


is valuable as a test of the principle of conservation of mass, as well 


as an independent check on the accuracy of the finite difference 


solution. 
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Further information on reservoir simulator development can be found 
in publications by Craft and Hawkins, 1959; Peaceman, 1969; Ford, 1971; 


Thomas, 1982; and Farouq Ali, 1980. 


Computational Procedure 

A computer program to solve the numerical model for a petroleum 
reservoir, as developed in the preceding section, was written and coded 
in FORTRAN IV language, and run on an Amdahl 470 computer under the MTS 
operating system. The program is composed of a series of short, well- 
defined sections, each responsible for solving a specific part of the 
problem, and contains two subroutines to allow for interpolation of data 
by cubic splines. The program may be briefly described as follows: 

Initially, the program reads in and prints out all required 
reservoir parameters, time-invariant data and program flags that assign 
various options to a particular simulation. Thereafter, the constant 
parts of the transmissibility terms are calculated (with the no-flow 
condition at the outer boundaries of the reservoir being represented by 
zero transmissibilities at the boundary points); the initial estimates 
of fluid properties are made; and pressure and saturation distributions 
are produced. The first major loop is then set up, and provides for 
advancement to the next time level when a solution for the present time 
has been obtained. The remainder of the transmissibility terms, as wel] 
as the coefficients Z, B, D, F, H, and S, are subsequently calculated. 
The second major loop is started to provide for further iteration when 
convergence to a solution has not been attained. Fluid properties are 


updated to the extrapolated time level, and the E coefficient and 
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qi j,k terms are determined. After solving for the pressure and 
saturation distributions, a check determines if the solutions have 
converged. If so, a mass balance check is performed, and if acceptable, 
the pressure and saturation distributions, together with the 
production/injection rates and mass balances, are printed out. A time 
increment is selected; all reservoir properties are updated; pressure 
and saturation distributions are extrapolated to the next time level; 
and the program returns to the outer loop. If the tests for convergence 
or the mass balance checks are unacceptable, the program updates 
properties to the latest estimates and returns to the inner loop. A 
maximum number of iterations is prescribed to account for the 
possibility of an oscillating solution, and if exceeded, the program 
breaks the loop to print a warning message. 


A flow chart of the computer program is presented in Figure 8. 
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Figure 8. Flow Diagram of Computational Procedure 
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Figure 8. (Continued) 


SHES 


\ 


a 


Ve 


etnsne!2 Spnapett: retelusiad 
anare ws seve 46 


— 7 


. Bile des 2aqtStit terug), Se 
| cing Pw an. Us cit i stind 


—am ese —_ 


ss. ap yrl sitrstec 
ay i alalaedl ) s 


- —————— 
Hire | ant a hee lila.) ; 
E akan y" Rec, h 

P ~ : : 


| — 
ee; f o a wi 
| Sea yn | 


337 


Update Coefficients, 
() Pressures, and Saturations 


No 
Check Check 
Pressure and Maximum Number 
Saturation Solution No of Iterations for 


for Convergence Convergence 


No 


Check 
Maximum Number 
of Iterations for 
Convergence 


Yes 


Check 
lacs 


Yes 


Stop 


Print Pressure and 
Saturation Distributions, 
Production/Injection Rates, 
Mass Balances 


End of 
Simulation Yes 


No 


Increment Time Level, 
Extrapolate New Pressure 
and Saturation Distributions 


Figure 8. (Continued) 
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CHAPTER FOUR 
LINEAR PROGRAMMING METHODS 
FOR AUTOMATIC HISTOR Y-MATCHING 
As noted in Chapter 2 (Review of the Literature), the first 
automatic history-matching algorithm that yielded feasible estimates of 
reservoir properties for practical field situations was developed by 
Coats et al. (1970), who employed a combination of least-squares and 
linear programming methods. While more complex automatic history- 
matching algorithms, often based on linear programming, were designed in 
later years, these did not significantly improve the original technique 
(Slater & Durrer, 1971; Thomas et al., 1972; Carter et al. 1974); and 
nonlinear optimization techniques could not guarantee attainment of a 
globally minimum matching error (Wasserman et al., 1975). Furthermore, 
although the system of partial differential equations, which constitutes 
the reservoir model is basically nonlinear, studies by Carter et al. 
(1970) support the contention of Coats et al. (1970) that the matching 
error-parameter relationship could be adequately approximated by a 
linear function, especially where reservoirs possess relatively high 
porosities and permeabilities (Slater & Durrer, 1971). 
Since the work of Coats et al. (1970) has attracted much attention 
- and most authors who considered it obtained supporting evidence and 
acknowledged the technique as at least a basically practical approach to 
automatic history-matching - it is pertinent to analyze the algorithm 


further in order to determine its validity more precisely. 
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40. 
Derivation of the Algorithm 
As developed by Coats et al. (1970), the least-squares, linear 
Programming algorithm for automatic history-matching is derived as 
follows: 
Kets x: 


J 
parameters that constitute the reservoir description data. In order to 


(j=1,2,...,0) represent the porosity and permeability 


reduce the number of these unknowns, the reservoir is partitioned into 
zones (based on existing geological data), which are each assumed to 
possess constant porosity and permeability. Now, let d; 
(i=1,2,...,1) represent performance data (pressure distributions, 


production/injection rates), with observed performance data denoted by 


anes and calculated data (from simulator runs employing various values 
of xj) denoted by eaa .- The matching error set, Ei > is defined 
by 

ane OPS ! qcale 


and the history-matching problem lies in determining a set of 
description parameters, % » which will minimize some norm of this 
error set. One notes here that, in order to utilize optimization 
techniques, the number of performance data must exceed the number of 
reservoir description parameters, i.e., [I > J. For each Xj to be 
estimated, it is also necessary to impose upper and lower bounds, 


denoted and Xj1> respectively (and based on existing geological 
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a 
A total of N simulations are then performed, each using a 
different set of random description parameters, Xj These random 


parameters, denoted xi for each run r (r=1,2,...,N), generate the 


qcalesr 


calculated performance data, ; 


» and the resulting matching 
r ‘ae : : 
errors, €, . For any run rr, each x is selected using a uniform 


random number generator, so that 


where R is a random number between O and lil. 

On the supposition that the matching error, Ei 5 is a single- 
valued function of the reservoir description parameters, Coats et al. 
(1970) chose to model this dependence by a linear relationship, and 
justified the approximation by pointing to the reasonable success of the 


resultant history-matching algorithm. They thus defined the relation 


between the matching errors and reservoir description parameters as 


’ (4.1) 


wheres 3)= 1,25... ,05, and | Xp ds 
A procedure for determining the unknown a,j; in Equation (4.1) was 
then developed by making use of a least-squares approach. For each 


A A le 5 : 
run sr, the deviation D. is defined as 
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N 
sum). ae is then minimized (differentiate and set to zero) to 


N 
2 [ xnxy Jay, = i ae , (4.2) 


where, 1021 2,...5). 0nd n= (0.1.2... Js, An alternative dernivation 
of Equation (4.2) is provided in the Appendix. 

For each i, Equation (4.2) represents J + 1 simultaneous linear 
equations in the J+ 1_ unknowns AjQ24jz]s-+esaGy - It is helpful to 


expand Equation (4.2), for any i, by writing 
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Note that the coefficients of the a terms are constant with respect 


ij 
to the index i, and that this method essentially decouples the equations 
for each e. from those of every other Ei - 

To solve for the unknown a;;, define a; and cj; as the column 


vectors (daasdyisese sds )). and (GieeSiis Cau respectively, and 
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the matrix B containing the Bhi coefficients, where 
i her roe 
a ee Ge ee ete en es 
nj eile! in Te 


a. = Bol. ; (4.3) 


anderepeaved applicationvor equation) (4.3) t0m1l=)l.2s.05s1 yields 


jc Since the Dn coefficients are independent 


of i, the matrix B need be calculated only once and stored. Note 


the entire set of a, 


that if the matrix B is found to be non-invertible, it is necessary to 
rerun the random number generator until an acceptable set of values is 
obtained. 


Equation (4.1) gives 


where the a are now known, and the inverse problem of determining a 


zig 
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best set of reservoir description paramters, Xj » 1s therefore reduced 

to minimizing some norm of the errors subject to the constraints 
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with the goal of eliminating negative porosity and permeability 


estimates. 
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44, 
In order to utilize a linear programming solution technique for 
minimizing Equation (4.4), CF is expressed in terms of slack 
variables: 


Sipe cpl cJeian cone! 
Combining this with Equation (4.1) then yields 
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wnere™ 1°=91,2,...,). Also; Equation (4.4) becomes 
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and by the triangle inequality, an upper bound on S_ is minimized when 


I 
Pehotal i 


I 

L, i 

attains a minimum value. Note that, since 
I J I 
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upper and lower bound for S can be computed. The linear programming 


problem can consequently be formulated as follows: 
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The total of 21 + 3J variables are defined as 


Reservoir description parameters: 


Slack variables for errors, oe 
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Slack variables for upper bound constraints on the Xj! 
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Slack variables for lower bound constraints on the Xj 
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Solution by the simplex method, which then furnishes a best set of 


estimates for the reservoir description parameters, is straightforward. 
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Computational Procedure 

A computer program for implementing the automatic history-matching 
algorithm of Coats et al. (1970) was written, coded in FORTRAN IV 
language, and run on an Amdahl 470 computer under the MTS operating 
system. As the algorithm required a petroleum reservoir simulator, the 
program was designed to be compatible with the simulator developed in 
Chapter 3. 

The program reads the initial reservoir data, including the 
observed performance data, dae ,» and the constraints on the 
description parameters to be estimated. A loop, enclosing the entire 


Simulator, then uses a random number generator to select the reservoir 


porosities and permeabilities, x , and calculates the performance 


qcalesr 


data, j 


,» and the errors, ey » for a specified number of runs. 
Using the description parameters, xt » and the errors a , the 

matrix 8B .and vector c; are then set up, and the aij coefficients 
are calculated. Thereafter, the coefficient matrix of the constraint 
equations for the linear programming problem is determined, and solution 
by the simplex method subsequently yields the best estimates of the 
reservoir description properties. 


A flow chart of the computer program is presented in Figure 9. 
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Figure: 9. Flow Diagram of Computational Procedure 
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Figure 9. 


Set Up Matrix B 
and Vector Cc; 
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Set up Matrix Containing 
Coefficients of Constraint 
Equations 


Solve the Linear 
Programming Problem 
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49. 
Discussion 

Although the experimental results reported by Coats et al. (1970), 
Slater & Durrer (1971), Thomas et al. (1972), and Carter et al. (1974) 
indicate that the automatic history-matching method proposed by Coats et 
al. (1970) produces reasonable estimates of reservoir description 
parameters, it is pertinent to examine the theoretical foundation of the 
algorithm. 

Without benefit of an experimental test of the algorithm, the 
assumed linear functional dependence of matching errors on reservoir 
description parameters, and the subsequent reasoning for deriving the 
automatic history-matching procedure, may, a priori, be expected to 
yield poor results. 


Consider, first, the method by which the coefficients are 


aaj 
determined. In least-squares calculations, the assumption of an 


approximate linear behaviour implies that 


Ee e eiaeet 


d x) 

where a 1S a constant, x iS 4 set of reservoir description 
parameters, and qcalc(x) are the calculated performance data based on 
the description parameters provided. For the complex system of nonlinear 
partial differential equations that describe the flow mechanics of a 
petroleum reservoir, such simple linear dependence would, intuitively, 
appear unlikely. That experimental testing of numerous reservoir 
simulators under various practical field situations has shown the 


approximation and resulting linear programming solution to be generally 


valid (Coats et al., 1970; 
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Slater & Durrer, 1971; Thomas et al., 1972; Carter et alegwl9s/4) could 
be attributed largely to the relative insensitivity of a reservoir 
Simulator - i.e., a fairly wide range of reservoir description 
parameters will yield similar calculated performance data (Coats, 1969). 
Furthermore, according to Slater & Durrer (1971), a linear functional 
dependence appears to be particularly suitable when the reservoir 
possesses relatively high porosities and permeabilities. Note that, 
since the reservoir description parameters for each run, xt » are 


multiplied in pairs (see Equation (4.2)), it appears that there is some 


cancellation of the error in solving for the coefficients. 


a 
Secondly, consider the statement of the linear functional 
relationship between the matching error (for any point of observation) 


and the reservoir description parameters as given in Equation (4.1): 


where Xj (j=1,2,...,0J) are the reservoir description parameters, the 
matching error, ec, iS given by e = qoPs _ ycale , and d2>S and acale¢ 
are, respectively, the observed and calculated performance data at any 
point (which depend on the reservoir description parameters). 

For convenience, let x7 and x represent the column vectors 


containing the true and estimated reservoir description parameters, 


respectively. Then the model on which Coats et al. base their algorithm 
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where a is the column vector of coefficients to be determined, 6 is 
a random error, t is the independent time variable, and Yeerepresents 
all other independent and dependent reservoir parameters not included in 


Xt and x. Rewriting Equation (4.5) yields 


obs ] 
( Phe a 


do (Xpoyst) = Xyo¥ot) + a'xp + 6. (4.6) 


Now, since an estimate of the reservoir description parameters, x, 
must be used to replace the unknown true parameters, x7, the model given 


by Equations (4.5) becomes, in practice, 


Ger ——, 


Xypo¥ot) - d Mast =a (4.7) 


which, ideally, should be nearly zero. Substituting Equation (4.6) into 


(4.7) yields 


Xsyst) + a'(xX7-x) +6 ~ (4.8) 
It is evident that this expression will, in general, and regardless of 
the value of the aij coefficients, approach zero only when x * xr. 
Because the algorithm of Coats et al. uses a reservoir simulator 
with several different sets of estimated description parameters, it is, 
of course, quite possible that the expression given by (4.8) yields a 


result close to zero. If the x’ are selected such that they are 
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ayer 
centred about the true reservoir description parameters, Xy, the error 
in estimating x7 can be balanced. Thus, if the bounds on the x are 
symmetrical about i.e. Ps eens ar ae ee 
y Xr (1.e., for each j, Xju 7 Xf = %; Xj1)> the 
random number generator will ensure a fairly even distribution of the 


x" about Xt» and the model presented in Equation (4.5) should be 
reasonable. 

Experimental results support this argument. As noted above, the 
success of the method developed by Coats et al. is well-documented in 
numerous case studies; and these indicate that the algorithm also yields 
satisfactory results for multiphase flow as well as when alternative 
zonation patterns are employed. 

In practical situations, it is not likely that bounds on the 
reservoir description parameters (estimated from available geological 
data) will be absolutely (or even approximately) symmetrical about the 
true description parameters. In that case, the algorithm of Coats et 
al. tends almost exclusively to choose as reservoir description 
estimates the extremal, bounding values. Results obtained by Dr. S.M. 
Farouq Ali and students at Pennsylvania State University (1973) evidence 
this characteristic problem. Working with a hypothetical, two- 
dimensional gas reservoir consisting of six Zones, it was desired to 
determine the absolute permeabilities. The available performance data 
and initial reservoir data, as well as bounds on the unknown 
permeabilities, are listed in Table 1, and the zonation is shown in 
Figure 10. (Note the asymmetry of the imposed bounds about the true 
reservoir permeabilities.) Results of five different applications of 


the automatic history-matching algorithm, and the true reservoir 


permeabilities, are provided in Table 2. 
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TABLE 1 

RESERVOIR DATA AND PERMEABILITY BOUNDS 
Welly Eocation: (2,2) 
Production Rate: 0.3277 m3/s 
AAT bilecks of equal size: Ax = Ay = 267 m 
Thickness: 9.14 m 
Temperature: 360°K 
Porositys 0.20 
Initial Pressure: 3447.4 kPa 


Time (days) Pressure at Well (kPa) 
100 SY. Aerw) 
200 3062.0 
300 2898.6 
400 2734.5 
500 25097 
600 2402.8 
700 220369 
800 2063.0 
900 1889.9 

1000 1714.0 
1100 1534.1 
Permeability Bounds Cum®) 


0.049 < ky < 0.295 


0.987 < ky < 4.935 
0.099 < kz < 0.493 
Delg7acmikne< 0200 
Deg evmcake 094022 


0.010 < kg < 9.869 
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Reservoir Zonation 
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Zone 
Zone 
Zone 
Zone 
Zone 


Zone 


Run 


6 


#1 


04296* 
2.061 

0.493* 
OAS iit 
05987* 


3.378 


ESTIMATED PERMEABILITIES (um 


#2 


0.049* 
1.433 

0.493* 
OMS 7= 
0.987* 
0.010* 


* indicates an extremal value 


TRUE PERMEABILITIES (um 


Zone 1 
Zone 2 
Zone 3 
Zone 4 
Zone 5 


Zone 6 


TABLE 2 


#3 


0.049* 
O59387% 
030998 
Oolg7s 
C298 7a 
O30 10% 


0.099 
1.974 
03099 
0.148 
2.467 
0.0001 
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#4 


0.049* 
459 35% 
0.493* 
05296* 
05987% 
0.010* 
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#5 


0.049* 
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Further analysis of the problem of asymmetrical bounds was carried 
out using the computer programs developed in Chapter 3 and in earlier 
sections of Chapter 4. Considering a hypothetical, two-dimensional, 
three-phase petroleum reservoir consisting of six zones, four 
applications of the algorithm of Coats et al. were used to estimate the 
absolute permeabilities. The first used bounds symmetrical about the 
true reservoir permeability; the second employed some asymmetrical 
bounds; the third used asymmetrical bounds entirely; and the fourth used 
intervals not containing the true permeabilities. The available 
performance and initial reservoir data and bounds on the unknown 
permeabilties are summarized in Table 3, and the zonation is shown in 
Figure ll. Results of the study are set out in Table 4. One observes 
that, as stated earlier, the history-matching algorithm selects extreme 
values the more frequently the less symmetrical the bounds on the 
permeabilities become about the true reservoir permeabilities. 

Coats et al. noted that, even in the case of symmetrical bounds 
(which they and all others exclusively employed in tests), the algorithm 
often chooses at least one or two reservoir description parameters at 
their upper or lower limits. They therefore suggested that repeated 
passes could, if desired, be made with shifted limits on the necessary 
parameters until all estimated values lie within specified bounds. 
However, they point out that it is far more convenient to choose bounds 
representing a reasonable range, apply the algorithm once, and accept 


any resulting extremal values selected. 
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TABLE 3 
RESERVOIR DATA AND PERMEABILITY BOUNDS 
Relative Permeabilities and Capillary Bounds 


0il-Water 0i1-Gas 
Sw Krew = Krow Pcow(Pa) “0, Keg Krog = Pego (Pa) 
Oct 0 1 28344 Cheuk OE 5.20 1 31144 


Oe2y 30.0016..0.375 655 0.2 0.410 0.009 462 
Dire eU 0S lien 0 57.35 496 Dom BU Sh Oe aU OSM 290 
Gee nd.0c59 940).590 421 O-4 >) 052205 80.062 138 
Oe Sieh W06/2.05420 352 Gooey UStt0e ) OUs 110 -7 

O26. O51000. 0.210 283 0.6 0.080 0.190 -152 
OeAee0.1400 0.070 214 Cao 0030p 80355 -297 
O50) 052000" 05016 145 053°) 0.00525 0237.0 -44] 
Useo.2500- —0 76 One 9 ee 1 -586 


where k.,, is the relative permeability of oil in the oil-water system, 
and Krog is the relative permeability of oil in the oil-gas system; 
and Sg = ONS 0 


phase gas-oil system, with irreducible water present. 


- So » where S._ is the oil saturation in the two- 


Formation Volume Factors, Fluid Viscosities, Gas Solubility 


(kP.) B,(m3/sm°) B, (m3/sm>) Bg (m>/sm>) Uy (mPa ,s) Rew(sm>/sm>)Reo (sm>/sm°) 


p 
0 1.0000 1.0000 1.0000 3.700 0 0 
1379. 11160 1.0000 0.0845 3.050 2 13.53 
2758 =: 1.1245 0.0413 2.879 2 19.77 
AiG VeiS75 0.0270 2.622 2 25.11 
5516 1.1500 0.0198 2.470 2 30.10 
6853 ~—- 1.1623 0.99354 0.0156 2.390 2 34.91 
11032. 1.1569 0.0091 2.545 2 34.91 
15168 1.1569 0.0065 2.700 2 34.91 
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TABLE 3 (continued) 


Ug 0.0167 mPaes; We 0.43 mPaes 

So Wesed tie Tee 

Initial porosity: 0.09 

Initial pressure: 3000 kPa 

initial) Oil ssaturation: 07/0 

Initial water saturation: 0.30 

Py sc = 806 kg/m? ; Byaee = 110g kg/m> ; Pee rere kg/m? 
Wet ocatl Ons l.4 smo.) uo sc) t Coaons (562 elo. 0 a Osa) 
Oil production rate in each well: 9,2x107° m/s 
THiEKNESS<genlcam 

Depth of block (1,4): 2060 m 

Dip: 7° downward in positive x-direction 

Bubble point pressure: 6853 kPa 


Dimensions (m) 


Ax, = 1500 Ay, = 500 
AX = 1200 Ay5 = 200 
AX3 = 2000 INE 500 
Axqg = 500 Ayq = 100 
Axs = 500 Ays = 1500 
Axe = 200 Aye = 100 


Performance Data: 272 days; N= 12 runs 


Well Block Pressure (kPa) 
(1,4) 2822 
col 2975 
(3,3) 2991 
(3,5) 3007 
(5,2) 3169 
(5,6) 3133 
(6,4) 3408 
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TABLE 3 (Continued) 


Permeability Bounds (un?) 


Case 1 (symmetrical) Case 2 (partially symmetrical ) 
Oso. < k, < 0.9 ORS ee ky 4 02) 
OAS ko < 1.0 0.7 < ky < 1.2 
O.1 < k3 <-057 0.1 < k3< Og7 
ers ies 15 Sse ee 
Oe ee, Oe 2 ee ee 
O.1 < ke < 0.9 0.4 < ke < 0.9 

Case 3 (asymmetrical) Case 4 (non-inclusive) 
0.2 < ky < 0.7 0.7 < k) < 1.5 
0.7 < ky < 1.2 0.2 < ky < 0.6 

0.35 < k3 < Ov7 0.8 < kz < Le 0 
0.255 ky s Led 0.2 < kg < 0.8 
0.1 < kp < (We) 1.0 <k,; < 1.4 
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ESTIMATED AND TRUE PERMEABILITIES (um) 


0.620 
10F 

0.376 
0.921 
0.819 


0.483 


* indicates an extremal value 


TABLE 4 


0.4* 


61% 


True 
0.6 
0.8 
0.4 
0.9 
0.8 


0.5 


62. 

In practice, it seems therefore appropriate to use the automatic 
history-matching algorithm of Coats et al. (1970) for single-phase flow, 
and multiphase flow where relative permeability curves are well 
defined. If it results in selection of most of the extremal values, a 
readjustment of the bounding values should produce reasonable data on a 
supplementary pass. The relative insensitivity of a reservoir simulator 
to small changes in permeability appears to enhance the acceptability of 
the results obtained by use of this algorithm. The reason for this is 
that small variations in parameter estimates will have little effect on 


future performance predictions made by the simulator. 
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CHARLER wei VE 
A DIRECT APPROACH TO AUTOMATIC HISTOR Y-MATCHING 

In 1972, Veatch & Thomas published an entirely new approach to 
automatic history-matching. While previous (and subsequent) methods for 
automatic history-matching resemble the algorithm of Coats et al. 
(1970), in that they are essentially ex post facto techniques (i.e., 
after numerous simulations using different initial data, various 
optimization and/or statistical methods based on comparisons of the 
known performance history with the results of the simulations are 
employed to determine acceptable estimates of the reservoir description 
parameters), Veatch & Thomas developed an easily implemented direct 
method for inverse reservoir simulation that substantially reduces the 
demand for computer time. 

By treating the finite difference analogues of the partial 
differential equations which model multiphase flow as a system of 
equations in the unknown reservoir description parameters, and employing 
performance histories over several periods of time as input data, it is 
possible to determine directly the description parameters in a single 
computer run. The particular advantage of this method is that it does 
not require implementation of a reservoir simulator, and thus eliminates 
computational inaccuracies arising from inherent inadequacies of the 
simulator. It should also be observed that the equations are solved 
directly, and that the solution method is therefore applicable to 
multiphase, compressible flow in heterogeneous reservoirs. Additionally, 
Veatch & Thomas included provisions for constructing ‘best' estimates of 
reservoir description parameters in cases where no unique solutions can 


be computed. 
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64. 
Veatch & Thomas' direct method for automatic history-matching 
represents a considerable saving in time and effort in handling such 
problems, and appears to have the potential to produce superior 
estimates of reservoir description parameters. Since, as indicated in 
the literature review, this procedure has so far attracted little 
comment, it is a matter of some interest to analyze it here as an 


alternative to Coats et al.'s algorithm. 


Derivation of the Algorithm 

Using the standard technique for representing a petroleum reservoir 
as a grid-block system, any partial differential phase equation required 
for describing the flow mechanics in any block (i,j) in a two- 


dimensional areal reservoir may be written in the form 


n n n 
OK CUT oe] mens, ehh pos SUT eM spp eaie OR Aer tire 
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N and 


where the permeabilities and porosity are unknown, and the ey 
P" are coefficients containing relative permeabilities, formation 
volume factors, fluid viscosities, spatial increments, and potentials, 
all at time level n (cf., for the case of three-phase flow, the 
equations derived in Chapter 3.) These coefficients, as well as the 
q’ terms, depend on both the pressure and saturation distributions 
throughout the grid-block system at several points in time; and since 


the only firm sources of information lie in blocks containing production 


or injection wells, various interpolation schemes must be employed to 
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Oo. 
estimate the required pressure and saturation distributions in the 
remainder of the reservoir. 

Theoretically, a solution to any one phase equation should satisfy 
all other phase equations, and Veatch & Thomas consequently based their 
algorithm on a single, arbitrarily chosen phase equation. The use of 
Equation (5.1) for any phase, when combined with known reservoir 
pressure and saturation distributions at several points in time (which 
must be at least as numerous as one plus the number of unknowns in each 
block) thus results in a linear system of equations in the unknown 
reservoir description parameters. 

Veatch & Thomas offered a detailed discussion of the matrix form of 
the resulting linear system of equations (which is of the familiar Ax = 
b form), and suggested solving it in a least-squares sense by applying 
Householder transformations (see Householder, 1965). They observed that 
Singular, nearly-singular, and ill-conditioned systems of equations led 
to non-unique solutions, and identified such cases from analysis of the 
diagonal elements of the matrix of coefficients. 

Moreover, Veatch & Thomas chose to assign permeability values to 
the faces, rather than the centres, of each grid block (i.e., ky was to 
be defined at (7+1/2,j) and (i-1/2,j) rathersthan at (1.3); ‘and’ the 
history-matching problem consequently contained five, rather than two, 
unknowns in any block. Given the large number of unknowns over an 
entire reservoir grid system, they suggested that only small groups of 
blocks be treated initially in order to construct local solutions. And 
since the most accurate pressure and saturation information on the 
reservoir lies in blocks containing production or injection wells, they 


further recommended that the initial computations be focussed about the 


well blocks. 


+o seu aHT nov Tune. semua: NSRON \Ahiet? tate vafpnle a Go init Progts 
stoviwagn word aotw ben tdmos ned eesiq yma nom (5.8) not tenes, 
estdw) enry nt gantoy (saver Je enortuctzerh notaeniiee bh emveeeTg 
dose at 2nwondne Ww redmue sd eT Sto en BUOTAMN ae Teeeh teed Seber 
awnieiny eno ni anotinups To mayaye steail & Af 2a tues zona (aooTd 

. star nos Sqtoseat aioyv1se8" 

io weet? Ktvenn Ot 10 Molegnoeth bel iecsh nubs Me aenonl ee oaeee 6 
© Kh sei(tan? ete Toet dain) abotraups, teomsreta Naantt ign a iveey GAS 
antyfqus ed sense ecisipe-dzeo7 Bat 3! patwine besespoue ‘baw (wet @ 
rot Hewrecho ven! . (20 ,sebfotladuoh dee} zn feeb tete ss sabtodgsauon 
bal encitaine 1® emageve benchtPbAose Tih ons reliwn! ee ineem letapate 
ait Yo 2teylens mo1) 28262 fouz bariisrabt bis .2nofsvine Saprnienan og 
emnsto}t*soa Yo xiatan arts Yo 2tnsmefe fenopeth 

of? aoulny v7! (idésened nefess Oo sane enmon? @ totes vevosom 

oo 2aw yd eet) Sole bino dose Foy 2astne oi? nid? vemeen _eecn? rtd 
aft tne fi(ek) fener Setar (iy Sidst) bre (Sieh) te bent tep ad 
Land fie sation .sver bseteitos yMtipsupencs ovidomm anédadamynodetht 


6 FaVO zONUNNNY Yo TSdMON Spe SL ovr .AvOld ys NE emeontAY 


ty equose Pen2 ino ton) beveeopue yet? ,wereye bine vtowaee ortgne 
stngtivfoa Teoo) tnuisenod GF Ihr Nt El TSrTAT betes s¢ 2k30fd 
sAt m NOtsaMIOs notputer BAG siue2ory atewoon. teow ed eonte 
wand )T Few NotszeEni ono zoubINd BART! asNy eield a ” Vovrsee7 
St’ 2und6, tise anv0" ad 2noPSeTumMO> fstitot ot? IeAt 


= 
= 


66. 
Thus, Veatch & Thomas' history-matching algorithm begins by making an 


initial pass over the entire reservoir in order to establish the systems 
of equations about each well block. Each system of equations is then 
tested for uniqueness; and if a unique solution can be computed, the 
reservoir description parameters are calculated. The systems of 
equations associated with the periphery of each block for which a unique 
solution was found are thereafter determined. If these systems produce 
unique solutions, the reservoir description parameters are computed, and 
the areas are gradually and uniformly enlarged over the reservoir until 
all possible unique solutions are obtained. If no unique solution can 
be found, the algorithm begins a non-unique construction of the 
solution. The computations in that case resemble the previous automatic 
history-matching algorithms, in that they require the imposition of 
upper and lower bounds on the reservoir description parameters, and are 
based on least-squares techniques in which the error in estimation is to 
some degree minimized. 

In tests of their algorithm on several hypothetical reservoirs, 
Veatch & Thomas found generally excellent agreement between actual and 
computed reservoir description parameters in the case of incompressible 
flow, for which unique solutions could be obtained. But when non-unique 
solutions were encountered, they conceded that "there is no assurance 
that the results of the prediction phase will be reliable", and, “as the 
prediction is made farther into the future, the results become even 
worse." The one example of compressible flow treated by Veatch & Thomas 
led to a non-unique solution, and implementation of a bounded least- 
squares technique produced reasonable results. However, as in the case 


of the automatic history-matching algorithm of Coats et al. (1970) the 
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apparent success of the history-match was due to imposition of 
symmetrical bounds about the true values of the reservoir description 


parameters. 


Discussion 

While Veatch & Thomas' method for automatic history-matching is 
potentially useful, it contains several inherent weaknesses. Its main 
disadvantage is its tendency to be very sensitive to the input pressure 
distributions. (Such sensitivity to minor variations in reservoir 
pressure was observed as early as 1964 by Jacquard.) This is 
illustrated by a simplified version of the direct method applied to a 
hypothetical reservoir. 

Consider the oil phase equation appropriate for modelling three- 


phase flow (Equation (3.7a)): 
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(Sac) 
For convenience, suppose the reservoir is two-dimensional, and 
Ges ky = k throughout the reservoir. Then, for any given block 
(i,j), Equation (5.2) may be written by finite difference techniques as 
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-VWallog) -(g) }-ag . (5.3) 


where the various terms are as defined in Chapter 3. 

To simplify the problem further, assume that the porosity is known, 
so that Equation (5.3) contains the permeability as the only unknown. 
Knowledge of pressure and saturation distributions at two points in time 
will then suffice for estimating the permeability in each grid-block 
throughout the reservoir. 

A three-phase, two-dimensional areal reservoir, divided into a 7x7 
grid-block system, was employed in the sample calculations. All 
necessary input data, including the true reservoir permeabilities and 
the pressure and saturation distributions at two time points, are shown 
in Table 5. To test the sensitivity of Equation (5.3) to the pressure 
distribution, the pressures at the first time step were perturbed in 
various blocks by increments as small as 6.895 kPa (1 psi). When the 
pressure values in each block (i,j), i = 1,2,3; j = 1,2,3,4 were 
increased by 6.895 kPa, the resulting permeability estimates in those 
and surrounding blocks were strongly affected, with some estimates 
displaying shifts in the range of 0./7 wm? aay aes! ume, In fact, the 
solution generated several negative estimates of permeability. 

It is possible to provide an heuristic argument to demonstrate the 
sensitivity of Equation (5.3) to pressure distributions. Consider the 


coefficient on the left hand side of Equation (5.3): 
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TABLE 5 
RESERVOIR DATA 


Helm OCac ions (lei yem(4 cet 4 4 | (oer 

Oil production rates: 9.2x107° m3/s inebiocks (lem); (454087 57) 
1.91077 m3/s in block (4,4) 

Exeluded blocks 1,1)e4721) 

All. blocks offeqtal sizesi Ax =Paye= 3025%m 

Thickness: 12.2 m 

Poroswty: 0.2 

Depth below sea level: 305 m 

Initial oil pressure: 13100 kPa 

Permeability: 0.296 am? 

Initial o1l-saturation: 0.550 

[Initial gas saturation: #05050 


Oil Pressure Distribution (60 days) 
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Oil Saturation Distribution (60 days) 
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TABLE 5 (Continued) 


Water Saturation Distribution (60 Days) 
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Pressure Distribution (100 days) 
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Oi1 Saturation Distribution (100 days) 
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In those parts of the reservoir in which the pressure distribution tends 
to be uniform, this expression will approximate zero. Since each Kaj 
is determined by dividing this value into the right hand side constant 
of Equation (5.3), it is evident that small changes in the value could 
have a significant effect on the outcome of the history-match. This 
point was noted by Veatch & Thomas: "The greatest deviations occurred 
in the farthest regions from the wells where the pressure distributions 
from one time to the next are not significantly perturbed. As a 
consequence we begin to see the effects of ill-conditioning where small 
errors in the coefficient matrix are greatly magnified in the solution 
vector." 

The sensitivity of this direct history-matching method to the 
pressure distributions was further corroborated during personal 
communication with Dr. G.W. Thomas. 

Although the direct method for automatic history-matching is 
advantageous in that it performs independently of a reservoir simulator, 
its sensitivity appears to render it ineffective as a viable approach to 
automatic history-matching in practical situations, where utilization of 


various interpolation schemes can lead to widely different estimated 


pressure distributions. 
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CHAPTER SIX 
SOME MATHEMATICAL CONSIDERATIONS 
Theoretical Suitability of Darcy's Equation of Flow 

As discussed in Chapter 3, mathematical models simulating 
multiphase, multi-dimensional fluid flow in a petroleum reservoir are 
ultimately based on Darcy's equation of flow (used in conjunction with 
the law of conservation of mass). The effectiveness of this equation as 
the fundamental relationship governing fluid flow behaviour in porous 
media has been established by experimental testing, and a large body of 
Supporting material (applied to actual field situations) is readily 
available (see, for example, Muskat, 1949; Collins, 1961; Scheidegger, 
1974; or Thomas, 1982). But equally important for the significance of 
the mathematical model is that Darcy's equation of flow has some 
theoretical foundations. 

In 1975, Houpeurt described a method - first published in 1955 in 
Revue de ]' Institut Francais du Pétrole - by which Darcy's equation of 
flow could be deduced from theoretical considerations of the Navier- 
Stokes equations of motion. 

In Cartesian coordinates, the Navier-Stokes equations of motion for 
an incompressible fluid with constant viscosity may be written as 
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where u, v, and w are the velocities in the 


directions, respectively, F is a vector of body force density, and 
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The law of conservation of mass (or the equation of continuity) is, for 


this case, 


du 
Ox 


For steady-state fluid flow, 
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Furthermore, assuming flow to be isothermal, laminar, viscous, and one- 


dimensional along the x-axis, the equation of continuity in the absence 


of external body forces is simply 
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and the Navier-Stokes equations reduce to 
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dx ay’ 37° 


(6.0) 


Houpeurt then dealt with two particular cases of flow. In the 
first instance, he considered flow along the x-direction through a three 
dimensional semi-infinite strip, defined by O0< x<L,-~¢ y¢o, 


-e/2 < z <e/2 , so that Equation (6.1) could be written as 


2 
d d 
cee ae (6.2) 
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dzé 
Since the two terms in Equation (6.2) depend upon different variables, 


each term may be set equal to a constant and integrated. For the 


boundary conditions 
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and 


Calculation of the output of fluid per unit time, q, over a distance 


b along the strip yields 


7.e., this is the fluid output per unit time through the cross-sectional 
area be. For a porous medium with porosity o , the total output of 


fluid per unit time across the strip, Q , is then given by 
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where A is the total cross-sectional area along the strip. 
In the second instance, Houpeurt considered flow through a 


cylindrical porous medium, and, by methods similar to those used above 
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where r is the radius of the cylinder. 

Comparison of Equations (6.3) and (6.4) then led Houpeurt to 
Suggest that, for any cross-sectional area A of a porous medium 
through which the direction of flow is normal, the effect of a pressure 


gradient on the output of fluid per unit time, q, would be of the form 
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where k is a constant. This is, in fact, the result obtained 
experimentally by Darcy in 1856. 

The preceding arguments allow the observation that the problem of 
obtaining a fundamental expression for fluid flow in a porous medium 
reduces, essentially, to determining the macroscopic behaviour of a 
porous medium which exhibits microscopic heterogeneity. In other words, 
since the problem lies in the existence of two widely differing length 
scales, or ‘scales of variation' (the smaller relating to heterogeneity, 
and the larger to macroscopic behaviour), the objective must be a 
characterization of the large scale behaviour of the medium by 
eliminating small scale variation from the known equations which 
describe flow behaviour on a small scale (i.e., the Navier-Stokes 
equations of motion). 

In recent years, several systematic, rigorous methods for dealing 
with such problems have been developed. These are based on various 


averaging techniques, including spatial, temporal, and stochastic 


ata ¥o wri s 
oF “Seaquol pst nade’ (beB) Gnb (£8) eno} 75up3: Yo. maatysqmed 
mutbsm 2yovsg 6 ‘jo. A éers Tenétss92-z2079 ye OT tot seoppue 

avyreeig & % 3a9Ve saz, femmon 2f wolT To norsaavth ote fiotiw dgoowt 
me} qt tes hfvow co SaneY aTnu seq btu to. sughiowdd ne aewthe 


Stir. 


tentudde rillderfanh ae) nb hi Sehr? LagavEnes gah 4 Stele 
eeatent Vaae ai teers reqnee 

Te welding Sie 2ed7 HOtIevrstde ant wotts atieeunns gatpaget 6aT 
net pal 2o7hq 8 AE do bit “S01 dotzeprets Ietnanebaut 2 gniaietda 


a 


5 To tui vEnss atabsenrger of? anintmetas ot .yitataneces . eee 

.tievow wade ol wad teeepoverer Sigodennatr rtididen doldw mule evorvog 

ntgnst gabrstteh ybabty on 30 eines 4 Jud nl 26f! msfaerg aay eonte 

. \¢! heragiisgar co pitisten 4s Agl {we em an iokivey ts eabene’ 4e .eetepe 

| § 40 2240 Sy trasude ony tute fod Sfawtoreae ut rope! off tne 7 
ye. om’ben SAF to sno velle? Slese pis! ant norsextoasninetia. 

dali enhT Tends nwonw 2b WONT dT t6M 5. sieoe Tem potsantatta| 

pe se-1B?veN sty ).9.1) ae tees bs ® svotveted wel) adirpesb 


-{norion to ee ttsups 


lersver , 2.88% IapoeY x” 


Pi. 
averaging. In 1977, Keller developed the so-called 'two-space method' , 
or ‘method of multiple scales', as means for deriving simplified 
equations. The method can handle large amplitude variations in the 
coefficients which may occur over a small length scale, and Keller 
clearly illustrated the technique by applying it to a boundary value 
problem for a second order partial differential equation. 

In a 1980 paper, Keller further illustrated the use of the two- 
Space method for the simple case of heat conduction in one dimension, 
and ohn used this approach to derive Darcy's equation of flow from the 
Navier-Stokes equations. The following is a brief description of the 
two-space method as applied by Keller to flow through a porous medium. 

Consider an equation with coefficients that vary on a small spatial 
scale. These coefficients are first represented by writing them as 
functions of x/e , where x is a position vector and « a small 
parameter proportional to the length scale. That this will cause a 
function f(x/e) to be subject to rapid variation is evident by noting 


that its derivative with respect to x, i.e. 


is large for small values of e« , even though f' may be bounded. 

For a flow of a compressible, viscous fluid through a porous 
medium, let e denote the ratio of the pore diameter to the macroscopic 
scale. The relationships governing flow behaviour on the microscopic 
SCAlCe 1 e@cnm nates dn teml Osta tne pores) are then the Navier-Stokes 


equations of motion, the equation of continuity, and the equation of 


state, namely: 
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o(S + ueVu) = -Vp + u(v + 307+ )u te fh (6.5) 


’ (656) 


and 


pts ‘p( Diem (6.7) 


where op, p, U, and wu are, respectively, the fluid density, pressure, 

velocity, and viscosity coefficient, and f is the external body force 

per unit volume. In addition, it is assumed that u = 0 on the surface 
of the pores, to account for viscosity effects. 


Keller first introduced the variables y and +t , defined by 


Wo = eens 


t= t/el/? a 


and wrote u, uy, p, p, and f in the forms 
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and 


f= f(xXey ete) 
Note that yw is small of order ge » 1.e., the viscosity decreases 
proportionally to the decrease in pore size. Also, the powers of ec 
employed vary among application of the two-space methods to different 
problems. There is no formula to determine appropriate powers, and 
intuition or trial-and-error must be used to choose powers that provide 
acceptable final results. (For further references to this problem, see, 
@.g., Bensoussan et al. (1978), Chap. 2; Kogelman & Keller (1973); or 
Larsen (1975).) Replacing V_ by rans “1, , and substituting the above 


variables into Equations (6.5)-(6.7) yields 


+ BL(v+ev, )° + Hv,+e,)(Wytev,) «Ji <f AS (6.8) 


and 


0 = e(p) - (6.10) 


Assuming that U, p, p, and # have a regular dependence upon e , 


they can be expanded into 


Pin 2 
teense sst20 2% aisessadts MS ‘atre-to ha 4 mi Mo 
3) We z1dwou ahd beth cS84e shan i aeesraen, 98 sone 
thonettrh nt zbosrem, ad Ye ieee ak 8 


bik: Sesion sontraorgGs enbmn@iat Os-eTumoy of\et lenstl 
soiverg sors 2vawon sefions i ibeu Sd 7eum Acrng-bansTaae vo. otstegMe! Sg 
aoe cmetiony. etd te sesdetetior with TON) © ert Lins akan 
mo 1 (ETEE) Agri # somtzcO% Fh VER CHFET)- te ay emeevORnet QAgve 7 
wrovis. aed gvtdu ri dadne) de « eee! Ne, oo oathelqar QCRCeLy nected 
abtaty(t.a}-{e-4) — 


=) 
> © 


ere oe « fityranered + Bie 
(aa) | 1 Att yrange sal pikes + fen) js + 


(Rd) - he Tare bar od . Ea 


(or. ; We - ‘ 


, 4 ey. sonebgan ci uoen: ava Fw it 3. at wim 
cant athens wt 


pire 


wa: 


U(X,y¥,T,€) ea 


eee pat eon Olle) ar, 


Pea Pat SED y eta Ol a), 


and 


Cae eU,(Xs¥.t) + O(e) , 


(Goi) 


(6212) 


(6.13) 
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Substituting (6.11)-(6.14) into Equations (6.8)-(6.10), and equating 


coefficients of the lowest power of « in each equation yields 
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Similarly, equating coefficients for the next lowest power of 
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From here, Keller proceeded to determine the variables upon which 
each function depends, and after averaging (i.e., integrating over a 
large domain, D, of the fluid, dividing the integral by the volume V 
of D, and letting D and V tend to infinity), obtained two relations 


which generalize Darcy's law for nonlinear, time-dependent, compressible 


flows: 
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where the solutions to (6.15) and (6.16) are written as functionals of 


Py ands ot aes Spee in the form 
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and 


Py (xs¥5t) a B Xa Yoon at ool Da ° 


With these relationships in mind, Keller then considered several 
special cases which more closely resemble the usual form of Darcy's 
equation - flow of incompressible fluid, steady incompressible flow with 
constant viscosity, and steady compressible flow. In the case of 
Steady, incompressible flow, the simplified equation is the Darcy 


equation of flow in one of its common forms: 


where A(x) is the cross-sectional area normal to the direction of 
flow. 

Thus, based on the equations governing flow on the microscopic 
level, it is possible to obtain, from purely theoretical considerations, 


a characterization of large scale flow behaviour in a porous medium. 


Alternative Approaches to the Inverse Problem 

As shown in Chapter 2, petroleum engineers have devoted 
considerable efforts to the development of feasible automatic methods 
that solve the inverse problem. But analogous inverse problems also 
occur frequently in science, medicine, and engineering, and a large body 
of knowledge exists elsewhere in the literature (see, for example, 
Kagiwada, 1974; Carasso & Stone, 1975; or references provided by Sagar 
et al., 1975. More extensive reference listings may be found in Payne 


(1975) and Tikhonov & Arsenin (1977).) 
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A somewhat less complex version of the problem of parameter 
identification in an underground reservoir arises when considering 
groundwater flow. The equation modelling such flow is very similar to 
that which models single-phase flow in a petroleum reservoir; and 
numerous papers seeking to determine the spatially varying coefficients 
of permeabilty and porosity are to be found in journals of water 
resources and hydrology. It is interesting to note that most of these 
attempts follow direct approaches, rather than the ex post facto methods 
almost exclusively used by petroleum engineers. 

Proposed methods for solving the inverse problem (often for the 
Steady-state condition only) include: linear programming (Kleinecke, 
1971), automatic solutions (Emsellem & de Marsily, 1971), use of 
Subjective information (Lovell et al., 1972; Nutbrown, 1975), Galerkin 
solutions (Frind & Pinder, 1973), finite elements (Neuman, 1973), direct 
identification based on approximated derivatives of the dependent 
variable (Sagar et al., 1975; Yakowitz & Noren, 1976), and quadratic 
programming (Chang & Yeh, 1976). However, there is little discussion of 
the sensitivity of calculated solutions to the required input data and 
few practical case studies have been reported. 

A common feature of these direct algorithms is that they demand 
prior knowledge of the derivatives of the dependent variable (i.e., the 
pressure) over at least a portion of the reservoir. In two papers 
published in 1981, Richter studied the inverse problem for underground 
reservoirs using such input data, and analyzed possible solutions and 


solution sensitivity under various conditions. 
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In one paper, Richter (198la) considered the basic equation which 
models groundwater flow and flow of components in a petroleum reservoir, 


10%", 


where u represents the pressure, q is a source/sink term, a is a 
transmissibility term containing reservoir permeability (a@>0), B is a 
storage term containing reservoir porosity (8 > 0), and a and 8g are 
considered to be functions of the spatial coordinates. The inverse 
problem then centres on identification of a and Bg from observed 
values of u and q. 

Under steady-state conditions, the partial differential equation is 


reduced to the hyperbolic equation 
VaeVu + aAu = -Gg we Fe clerk” 


and in a systematic theoretical analysis of this problem, Richter 
employed a direct approach (i.e., a method involving an approximate 
solution of the hyperbolic equation). But for estimating a , this has 
severe practical limitations since it requires knowledge of the 
derivatives of u. 


Richter considered the problem under four different conditions: 


(i) inf ,|Vvu| UE as 


(17) inf Au > Oars 


ee ee hs «van aap “7 
avteyn| 287) addbeteshiun” tdcdage’ Ads, 1 Shotiaoneb let ot bensbtenas 
berreedo man} £ tnd, 9, Yo NOtss0/) /Aasbt) no epraas oott-maideng 

| 7 i 
et wot teens tstansnetttn, Tattinge aid ‘titi atinigiage ean 
Ret taupe ondommeqigt aff 03 bagube 


: ; a 
velath .mefdong' ata) fo Lesh Tertawiont? otdamndaye & at ibn 
svemrxougge Mm jn Fyn fa ey eS) nuvage tev! 6 beyolgne 
aphoahde we eabentiing 16° a) stool sties otoceqyh uta Ye nerzulor: 
afl to sobatane arivie sh a0 anottestntt Tagtwaeng. mpver 
7 Some Wo eeyigavtaed 
ihe hes: iayett in 907 so ae ae teneiteme namtatt 


$5. 


(iii) inf [max{|vu],au}] > 0, and 
(iv) inf [max {| vu] ,|Au| }] > dae 


Under each of the first three conditions, he showed that a unique 
solution a exists for any q , provided that the initial data 
(prescribed values along a portion of aQ ) are appropriately specified, 
and that a depends continuously on q , Vu, and Au. He also 
described the cases for which the fourth condition will lead to a unique 
solution. 

Examining the first condition, which implies the existence of non- 
intersecting characteristics, Richter was able to derive a bound which 
defines the sensitivity of a solution a to perturbations in q_ and 
ue He noted, however, that the bound is of little practical value, 
Since the measured pressures upon which the bound condition is based are 
not known with sufficient accuracy. As well, the bound suggests that an 
acceptable estimate of a is possible only if the observed u_ is 
sufficiently precise to allow accurate approximation of Au . 

The second condition allows for intersecting characteristics, which 
occur in the neighbourhood of a source or sink (a ‘point of 
degeneracy'). In this case, Richter showed that there is a preferred 
sense of direction along the characteristics, i.e., depending upon the 
initial starting point, the solution is drawn along the characteristics 
toward either the boundary or the point of degeneracy. 

Considering the third condition, which allows for the vanishing of 


the first order derivatives at points where Au is of the same sign 
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(thus, u cannot contain both maxima and minima in @ 1, Richter 
discussed the difficulty of prescribing q and boundary conditions on 
u which will guarantee that the forward solution satisfies |Vu| > 0 and 
Au > 0 throughout @. However, by combining results for the first two 
conditions, he derived a statement that provided for the possibility of 
a unique solution. It should be noted that some difficulty was 
encountered in reproducing the proof of this statement. 

Richter then proposed a particularly useful set of test conditions 
for measuring u , and, under these conditions, showed that a unique 
solution a exists for the hyperbolic problem without requiring Cauchy 
data. He also obtained a bound on the stability of the solution in 
terms of a,q =, and relevant properties of Q. 

In a companion paper, Richter (1981b) developed a finite difference 
method for approximating the solution a numerically, and proved that, 


under the third condition, the solution does in fact converge. 
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CHAPTER SEVEN 
SUMMARY AND CONCLUSIONS 

An intuitive, trial-and-error approach to history-matching can be 
costly and time-consuming. In an attempt to reduce these factors, 
considerable efforts have been made to automate the history-matching 
procedure for implementation on high-speed computers; and several 
methods based on ex post facto techniques have been developed for this 
purpose. 

This study examined presently available automatic history-matching 
algorithms reported in journals of petroleum engineering; discussed the 
viability of each; and explored the fundamental weaknesses of the two 
most promising methods. In addition, the existence of alternative 
approaches to solving the inverse problem for analogous, but less 
complex, situations was noted. 

As discussed in Chapter 2, the automatic history-matching 
algorithms available to petroleum engineers involve 

(a) iterative adjustment and regression analysis (Kruger, 1961; 

Jacquard & Jain, 1965; Jahns, 1966; Dupuy, 1968), 

(b) linear and nonlinear programming techniques (Coats et al. 
1970; Slater & Durrer, 1971; Thomas et al., 1972; Boberg et 
alee 978 iGanternety alatcil974 ja 

(c) energy dissipation analysis (Nelson, 1968), 

(d) optimal control theory (Chen et al., 1974; Chavent et al., 
1973; Wasserman et al., 1975; Watson et al., 1980), 

(e) direct solution (Veatch & Thomas, 1971), and 

(f) Bayesian estimation (Gavalas et al., 1976). 
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These publications generally treat only single-phase or pseudo 
Single-phase incompressible reservoir flow, although many would appear 
to be satisfactory for multiphase flow. A two-phase, incompressible flow 
problem was investigated by van den Bosch & Seinfeld (1977), but 
solution methods are suitable only for one-dimensional flow. Veatch & 
Thomas (1971) developed a direct method for automatic history-matching 
which is, theoretically, applicable to multiphase, compressible flow. 
However, for proper consideration of any multiphase flow problem, it is 
necessary to have knowledge of the relative permeabilities (which are 
utilized in the reservoir simulator of ex post facto techniques, or, in 
the case of Veatch & Thomas, directly in the phase equations being 
solved). While geological samples provide some rough estimate of these 
values, any history-matching problem should incorporate determination of 
the relative permeabilities as well as the absolute permeabilities and 
porosities. 

The most serious problem encountered in automatic history-matching 
is the tendency to construct ill-conditioned systems of equations (j.e., 
for the problem Ax = b, small relative changes in the matrix A or 
vector b produce large relative errors in the solution vector x). BY 
the very nature of the history-matching problem, inherent uncertainties 
exist in both A and b_ because they are based on a measured 
performance history. Further error inevitably arises from the finite 
word length of the computer. As Rust and Burrus (1972) noted: "the 
presence of these errors makes it impossible to obtain a meaningful 
solution by simply applying one of the classical methods to the system 


by itself. Moreover, if all we know is the system itself, there is no 
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nonclassical method that will give a meaningful solution" (see also 
Tikhonov and Arsenin, (1977)). The history-matching problem can 
therefore be treated only by incorporating some a priori information 
about the solution vector x. But due to the gross inaccuracy of any 
geological estimates (based on scarce physical data), even this 
requirement cannot be adequately satisfied. 

It is obvious, then, that current automatic history-matching 
methods are unsatisfactory. As expressed by Boberg et al. (1973), as 
well as by Dr. G.W. Thomas during personal communication, such methods 
are simply not competitive with manual techniques which an experienced 
reservoir engineer might employ; and practising reservoir engineers 
consequently continue history-matching on a trial-and-error basis, with 
good judgement and intuition providing reasonably reliable results. It 
would appear, then, that the determination of effective and reliable 
solutions to the inverse problem requires development of improved 
methods for estimating the geological structure and properties of the 
reservoir, and/or adaptation of sophisticated mathematical techniques to 


allow for a more thorough theoretical treatment of the entire problem. 
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APPENDIX 


Alternative Derivation of Equations (4.2) 


For any i, and a total of N runs, Equation (4.1) yields 


(1) (1) (a) 1) ly 
TOG dr ae ot ee 

(2) (2) (2) (2) pee) 
Tithe: Wa GUS Le @ be Fay eee a Ea oe 

(N) (N) (N) (N)-_ _{N) 
ean nee ©) 2 gee oe Ark, pees 


Now, multiply each of the N equations by its respective coefficient 
of asgs and add the resulting equations. Similarly, multiply each of 
the N equations by its respective coefficient of a.,, and add these 
resulting equations. Continuing this process for each coefficient of 
aij: j = 051,2,...,7, yields a total of wv + 1 Simultaneous equations 


identical to Equation (4.2). 
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